function solidTidesInfluence = SolidTidesInfluence(dateValues, APhi, ALambda, AHight)
  DayPerMonth = [31; 28; 31; 30; 31; 30; 31; 31; 30; 31; 30; 31];
  SecToRad = pi/(180*3600);
  GradToRad = pi/180;
  G_MMoon = 6.670E-8 * 7.3537E25;
  G_MSun = 6.670E-8 * 1.993E33;
  MaxIteration = 1000;
  
  APhi = APhi*pi/180;
  ALambda = ALambda*pi/180;
  
  for d=1:length(dateValues)
      vec = datevec(dateValues(d));
      datestr(dateValues(d),'yyyy-mm-dd HH:MM');
      AYear = vec(1);
      AMonth = vec(2);
      ADay = vec(3);
      AHour = vec(4);
      AMin = vec(5);
      ASec = vec(6);
      
      i=0;
      T1900=-0.5;
      for j=1900 : AYear-1
        if i==4 
            i=0;
        end
        T1900 = T1900+365;
        if i==0
           T1900 = T1900+1;
        end
        i = i+1;
      end
        
      for j=1 : AMonth-1
        T1900 = T1900+DayPerMonth(j);
      end
      
      if ((AMonth>=3) && (i==4))
        T1900 = T1900+1;
      end
      T1900 = (T1900+ADay-1)/36525;
        
      StarTime = 6*3600+41*60+50.548 + 8640184.813*(T1900-1);
      DayShare = 3600*AHour+AMin*60+ASec + (AYear - 1980)*0.8 + 50.5;
      StarTime = (StarTime + DayShare*1.002738)*2*pi/86400 + ALambda;
      T1900 = T1900+DayShare / 86400. / 36525.;
      T2000 = T1900-1;
    
      EclipticalSlop = (84381.448 -46.8150*T2000)*SecToRad;
      l = (296*3600+6*60+16.59 + 1717915856.79*T1900 + 33.09*T1900*T1900)*SecToRad;
      l1 = (358*3600+28*60+33.00 + 129596579.10*T1900 - 0.54*T1900*T1900)*SecToRad;
      F = (11.*3600+15*60+3.20 + 1739527290.54*T1900 - 11.56*T1900*T1900)*SecToRad;
      D = (350*3600+44*60+14.95 + 1602961611.18*T1900 - 5.17*T1900*T1900)*SecToRad;
      LGreat = (270*3600+26*60+2.99 + 1732564379.31*T1900 - 4.08*T1900*T1900)*SecToRad;
      
      Blunder = LGreat + (6.289*sin(l) - 1.274*sin(l-2*D) ...
        + 0.658*sin(2*D) + 0.214*sin(2*l) - 0.186*sin(l1) - 0.114*sin(2*F)...
        - 0.059*sin(2*l-2*D) - 0.057*sin(l+l1-2*D) + 0.053*sin(l+2*D)...
        - 0.046*sin(l1-2*D) + 0.041*sin(l-l1) - 0.035*sin(D)...
        - 0.030*sin(l+l1))*GradToRad;

      MoonLongitude = IntoCircle(Blunder);
      MoonLatitude = (5.128*sin(F) + 0.281*sin(l+F) - 0.278*sin(F-l)...
        - 0.173*sin(F-2*D) + 0.055*sin(F+2*D-1) - 0.046*sin(l+F-2*D)...
        + 0.033*sin(F+2*D))*GradToRad;

      MoonDistance = 384404.E5/(1 + 0.0545*cos(l) + 0.0100*cos(l-2*D) + 0.0082...
        * cos(2*D) + 0.0030*cos(2*l) + 0.0009*cos(l+2*D) + 0.0006*cos(l1-2*D)...
        + 0.0004*cos(l+l1-2*D) + 0.0003*cos(l-l1));
        
      [MoonLatitude, MoonLongitude] = EclipticalToEquatorial(MoonLatitude, MoonLongitude, EclipticalSlop);
      
      M = (1287099.804 + (99*360*3600 + 1292581.224)*T2000 - 0.577*T2000*T2000)*SecToRad;
      G = (1018578.046 + 6190.046*T2000 + 1.666*T2000*T2000)*SecToRad;
      e = 0.016730 - 0.0000419*(T2000+0.5);
      M = IntoCircle(M);
      Blunder = M;
      
      Old = 0;
      j = 1;
      while ((abs(Blunder-Old) > 1.0e-20) && (j < MaxIteration))
         Old = Blunder;
         Blunder = M+e*sin(Blunder);
        j = j+1;
      end

      SunDistance = 1.496E13*(1 - e*cos(Blunder));
      v = atan((sqrt(1 - e*e)*sin(Blunder)/(cos(Blunder) - e)));
      if abs(v-Blunder)>pi/8              
        v = v-pi;
      end
      SunLatitude = 0;
      SunLongitude = IntoCircle(v+G);
      [SunLatitude, SunLongitude] = EclipticalToEquatorial(SunLatitude, SunLongitude, EclipticalSlop);

      Radius = AHight*1E+2 + 6.378137E+8*(0.99832707 + 0.00167664*cos(2*APhi));

      cosZ = sin(APhi)*sin(SunLatitude) + cos(APhi)*cos(SunLatitude)*cos(StarTime - SunLongitude);
      Blunder = G_MSun/SunDistance*Radius/SunDistance*(3*cosZ*cosZ-1)/SunDistance;

      cosZ = sin(APhi)*sin(MoonLatitude) + cos(APhi)*cos(MoonLatitude)*cos(StarTime - MoonLongitude);
      Blunder = Blunder + G_MMoon*Radius/MoonDistance*(3*cosZ*cosZ-1+3*Radius*(5*cosZ*cosZ*cosZ-3*cosZ)/(2*MoonDistance))/(MoonDistance*MoonDistance);
    
      solidTidesInfluence(d) = - Blunder*1.2E6;
  end

function A = IntoCircle(A)

A = (A/(2*pi)-floor(A/(2*pi)))*2*pi;
% --------old version-------
%   if A>0
%       while A>pi
%           A = A-2*pi;
%       end
%   else
%       while A<-pi
%           A = A+2*pi;
%       end
%   end
%---------------------------


function[ALatitude, ALongitude] = EclipticalToEquatorial(ALatitude, ALongitude, AEclipticalSlop)

  s = cos(AEclipticalSlop)*sin(ALatitude) + sin(AEclipticalSlop)*cos(ALatitude)*sin(ALongitude);
  cc = cos(ALatitude)*cos(ALongitude);
  cs = -sin(ALatitude)*sin(AEclipticalSlop) + cos(ALatitude)*cos(AEclipticalSlop)*sin(ALongitude);
  atn = atan(sqrt((s*s)/(1-s*s)));
  if s<0
    ALatitude = -atn;
  else
    ALatitude = atn;
  end
  
  atn = atan(cs/cc);
  if cc<0
    ALongitude = atn+pi;
  else
    ALongitude = atn;
  end
